
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>odespy.problems &mdash; Odespy API 0.2 documentation</title>
    
    <link rel="stylesheet" href="../../_static/pyramid.css" type="text/css" />
    <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    '../../',
        VERSION:     '0.2',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="../../_static/jquery.js"></script>
    <script type="text/javascript" src="../../_static/underscore.js"></script>
    <script type="text/javascript" src="../../_static/doctools.js"></script>
    <script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
    <link rel="top" title="Odespy API 0.2 documentation" href="../../index.html" />
    <link rel="up" title="Module code" href="../index.html" />
<link rel="stylesheet" href="http://fonts.googleapis.com/css?family=Neuton&amp;subset=latin" type="text/css" media="screen" charset="utf-8" />
<link rel="stylesheet" href="http://fonts.googleapis.com/css?family=Nobile:regular,italic,bold,bolditalic&amp;subset=latin" type="text/css" media="screen" charset="utf-8" />
<!--[if lte IE 6]>
<link rel="stylesheet" href="../../_static/ie6.css" type="text/css" media="screen" charset="utf-8" />
<![endif]-->

  </head>
  <body>

    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="../../genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="../../py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="../../np-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li><a href="../../index.html">Odespy API 0.2 documentation</a> &raquo;</li>
          <li><a href="../index.html" accesskey="U">Module code</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <h1>Source code for odespy.problems</h1><div class="highlight"><pre>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">solvers</span> <span class="kn">import</span> <span class="n">compile_f77</span>

<div class="viewcode-block" id="Problem"><a class="viewcode-back" href="../../problems.html#odespy.problems.Problem">[docs]</a><span class="k">class</span> <span class="nc">Problem</span><span class="p">:</span>
    <span class="sd">&#39;&#39;&#39;</span>
<span class="sd">    The user derives a subclass and implements the right-hand side</span>
<span class="sd">    function::</span>

<span class="sd">        def f(self, u, t):</span>
<span class="sd">            &quot;&quot;&quot;Python function for the right-hand side of the ODE u&#39;=f(u,t).&quot;&quot;&quot;</span>

<span class="sd">    Optionally, the Jacobian can be computed::</span>

<span class="sd">        def jac(self, u, t):</span>
<span class="sd">            &quot;&quot;&quot;Python function for the Jacobian of the right-hand side: df/du.&quot;&quot;&quot;</span>
<span class="sd">    Some solvers also allow constraings (DAE problems)::</span>

<span class="sd">        def constraints(self, u, t):</span>
<span class="sd">            &quot;&quot;&quot;Python function for additional constraints: g(u,t)=0.&quot;&quot;&quot;</span>

<span class="sd">    Some problem classes will also define command-line arguments::</span>

<span class="sd">        def define_command_line_arguments(self, parser):</span>
<span class="sd">            &quot;&quot;&quot;</span>
<span class="sd">            Initialize an argparse object for reading command-line</span>
<span class="sd">            option-value pairs. `parser` is an ``argparse`` object.</span>
<span class="sd">            &quot;&quot;&quot;</span>

<span class="sd">    Other functions, for which there are default implementations</span>
<span class="sd">    in the superclass, are ``u_exact`` for returning the exact</span>
<span class="sd">    solution (if known), ``verify`` for checking that a numerical</span>
<span class="sd">    solution is correct (within a tolerance of the analytical solution),</span>

<span class="sd">    See the tutorial for examples on subclasses of ``Problem``.</span>
<span class="sd">    &#39;&#39;&#39;</span>
    <span class="n">stiff</span> <span class="o">=</span> <span class="bp">False</span>    <span class="c"># classification of the problem is stiff or not</span>
    <span class="n">complex_</span> <span class="o">=</span> <span class="bp">False</span> <span class="c"># True if f(u,t) is complex valued</span>
    <span class="n">not_suitable_solvers</span> <span class="o">=</span> <span class="p">[]</span>  <span class="c"># list solvers that should be be used</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&#39;&#39;</span>     <span class="c"># one-line problem description</span>

<div class="viewcode-block" id="Problem.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.Problem.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="k">pass</span>
</div>
<div class="viewcode-block" id="Problem.__contains__"><a class="viewcode-back" href="../../problems.html#odespy.problems.Problem.__contains__">[docs]</a>    <span class="k">def</span> <span class="nf">__contains__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">attr</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Return True if attr is a method in instance self.&quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="nb">hasattr</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">attr</span><span class="p">)</span> <span class="ow">and</span> <span class="nb">callable</span><span class="p">(</span><span class="nb">getattr</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span><span class="n">attr</span><span class="p">))</span>
</div>
<div class="viewcode-block" id="Problem.terminate"><a class="viewcode-back" href="../../problems.html#odespy.problems.Problem.terminate">[docs]</a>    <span class="k">def</span> <span class="nf">terminate</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">step_number</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Default terminate function, always returning False.&quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="bp">False</span>
</div>
<div class="viewcode-block" id="Problem.get_initial_condition"><a class="viewcode-back" href="../../problems.html#odespy.problems.Problem.get_initial_condition">[docs]</a>    <span class="k">def</span> <span class="nf">get_initial_condition</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Return vector of initial conditions.</span>
<span class="sd">        Not necessary to implement in subclass if the</span>
<span class="sd">        initial condition is stored in self.U0</span>
<span class="sd">        (this method then returns that condition).</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">if</span> <span class="nb">hasattr</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="s">&#39;U0&#39;</span><span class="p">):</span>
            <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="k">raise</span> <span class="ne">NotImplementedError</span><span class="p">(</span>
                <span class="s">&#39;class </span><span class="si">%s</span><span class="s"> must implement get_initial_condition&#39;</span> <span class="o">%</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">__class__</span><span class="o">.</span><span class="n">__name__</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="Problem.default_parameters"><a class="viewcode-back" href="../../problems.html#odespy.problems.Problem.default_parameters">[docs]</a>    <span class="k">def</span> <span class="nf">default_parameters</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Compute suitable time_points, atol/rtol, etc. for the</span>
<span class="sd">        particular problem. Useful for quick generation of test</span>
<span class="sd">        cases, demos, unit tests, etc.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="p">{}</span>
</div>
<div class="viewcode-block" id="Problem.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.Problem.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Implementation of the exact solution.</span>
<span class="sd">        Return None if no exact solution is available.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="bp">None</span>
</div>
<div class="viewcode-block" id="Problem.verify"><a class="viewcode-back" href="../../problems.html#odespy.problems.Problem.verify">[docs]</a>    <span class="k">def</span> <span class="nf">verify</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">atol</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">rtol</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Return True if u at time points t coincides with an exact</span>
<span class="sd">        solution within the prescribed tolerances. If one of the</span>
<span class="sd">        tolerances is None, return max computed error (infinity norm).</span>
<span class="sd">        Return None if the solution cannot be verified.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="n">u_e</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">u_exact</span><span class="p">(</span><span class="n">t</span><span class="p">)</span>
        <span class="k">if</span> <span class="n">u_e</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
            <span class="k">return</span> <span class="bp">None</span>

        <span class="k">if</span> <span class="n">atol</span> <span class="ow">is</span> <span class="bp">None</span> <span class="ow">or</span> <span class="n">rtol</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
            <span class="k">if</span> <span class="ow">not</span> <span class="p">(</span><span class="n">u</span><span class="o">.</span><span class="n">shape</span> <span class="o">==</span> <span class="n">u_e</span><span class="o">.</span><span class="n">shape</span><span class="p">):</span>
                <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;u has shape </span><span class="si">%s</span><span class="s"> and u_e has </span><span class="si">%s</span><span class="s">&#39;</span> <span class="o">%</span>
                                 <span class="p">(</span><span class="n">u</span><span class="o">.</span><span class="n">shape</span><span class="p">,</span> <span class="n">u_e</span><span class="o">.</span><span class="n">shape</span><span class="p">))</span>
            <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">u</span> <span class="o">-</span> <span class="n">u_e</span><span class="p">)</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">allclose</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">u_e</span><span class="p">,</span> <span class="n">rtol</span><span class="p">,</span> <span class="n">atol</span><span class="p">)</span>

    <span class="c"># subclasses may implement computation of derived</span>
    <span class="c"># quantities, e.g., energy(u, t) etc</span>
</div></div>
<div class="viewcode-block" id="convergence"><a class="viewcode-back" href="../../problems.html#odespy.problems.convergence">[docs]</a><span class="k">def</span> <span class="nf">convergence</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">u_ref</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">t_ref</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Given a series of solutions and corresponding t arrays in</span>
<span class="sd">    `u` and `t`, use the analytical solution or a reference solution</span>
<span class="sd">    in `u_ref` and `t_ref` to compute errors. Compute pairwise</span>
<span class="sd">    convergence rates for errors on consecutive time meshes.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">raise</span> <span class="ne">NotImplementedError</span>
</div>
<div class="viewcode-block" id="Linear1"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear1">[docs]</a><span class="k">class</span> <span class="nc">Linear1</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Linear solution of a nonlinear ODE,</span>
<span class="sd">    which is normally exactly reproduced by a numerical method.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;u&#39; = a + (u - a*t - b)**c, u(0)=b&quot;</span>

<div class="viewcode-block" id="Linear1.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear1.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">b</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="mi">2</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">c</span> <span class="o">=</span> <span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">c</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="n">b</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">not_suitable_solvers</span> <span class="o">=</span> <span class="p">[</span>
            <span class="s">&#39;AdaptiveResidual&#39;</span><span class="p">,</span> <span class="s">&#39;Lsodi&#39;</span><span class="p">,</span> <span class="s">&#39;Lsodis&#39;</span><span class="p">,</span> <span class="s">&#39;Lsoibt&#39;</span><span class="p">]</span>
        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">not_suitable_solvers</span> <span class="o">+=</span> <span class="p">[</span><span class="s">&#39;lsoda_scipy&#39;</span><span class="p">,</span> <span class="s">&#39;odefun_sympy&#39;</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="Linear1.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear1.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">+</span> <span class="p">(</span><span class="n">u</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">t</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">)</span><span class="o">**</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span>
</div>
<div class="viewcode-block" id="Linear1.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear1.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="o">*</span><span class="p">(</span><span class="n">u</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">t</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="Linear1.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear1.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">t</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span>
</div>
<div class="viewcode-block" id="Linear1.default_parameters"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear1.default_parameters">[docs]</a>    <span class="k">def</span> <span class="nf">default_parameters</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="k">return</span> <span class="p">{</span><span class="s">&#39;time_points&#39;</span><span class="p">:</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mf">2.</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="mi">3</span><span class="p">),</span>
                <span class="s">&#39;max_iter&#39;</span><span class="p">:</span> <span class="mi">10</span><span class="p">,</span> <span class="s">&#39;eps_iter&#39;</span><span class="p">:</span> <span class="mf">1E-6</span><span class="p">}</span>
</div></div>
<div class="viewcode-block" id="Linear2"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2">[docs]</a><span class="k">class</span> <span class="nc">Linear2</span><span class="p">(</span><span class="n">Linear1</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Linear solution of nonlinear 2x2 system,</span>
<span class="sd">    which is normally exactly reproduced by a numerical method.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;2x2 nonlinear system with linear solution&quot;</span>

<div class="viewcode-block" id="Linear2.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">b</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="mi">2</span><span class="p">):</span>
        <span class="n">Linear1</span><span class="o">.</span><span class="n">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">c</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="Linear2.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">+</span> <span class="p">(</span><span class="n">u</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">t</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">)</span><span class="o">**</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">,</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">+</span> <span class="p">(</span><span class="n">u</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">t</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">)</span><span class="o">**</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="Linear2.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="p">[[</span><span class="mi">0</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="o">*</span><span class="p">(</span><span class="n">u</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">t</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="o">-</span><span class="mi">1</span><span class="p">)],</span>
                <span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="o">*</span><span class="p">(</span><span class="n">u</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">t</span> <span class="o">-</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">)</span><span class="o">**</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="o">-</span><span class="mi">1</span><span class="p">),</span> <span class="mi">0</span><span class="p">]]</span>
</div>
<div class="viewcode-block" id="Linear2.spcrad"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2.spcrad">[docs]</a>    <span class="k">def</span> <span class="nf">spcrad</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">J</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">jac</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">)</span>
        <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">eigvals</span><span class="p">(</span><span class="n">J</span><span class="p">))</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
</div>
<div class="viewcode-block" id="Linear2.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">func</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">t</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span>
        <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">func</span><span class="p">,</span> <span class="n">func</span><span class="p">])</span><span class="o">.</span><span class="n">transpose</span><span class="p">()</span>

</div></div>
<div class="viewcode-block" id="Linear2t"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2t">[docs]</a><span class="k">class</span> <span class="nc">Linear2t</span><span class="p">(</span><span class="n">Linear2</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Linear solution of trivial 2x2 system,</span>
<span class="sd">    which is normally exactly reproduced by a numerical method.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;2x2 trivial system with linear solution&quot;</span>

<div class="viewcode-block" id="Linear2t.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2t.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span>
                <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="Linear2t.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2t.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="p">[[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">]]</span>
</div>
<div class="viewcode-block" id="Linear2t.spcrad"><a class="viewcode-back" href="../../problems.html#odespy.problems.Linear2t.spcrad">[docs]</a>    <span class="k">def</span> <span class="nf">spcrad</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="mf">0.0</span>


</div></div>
<div class="viewcode-block" id="Exponential"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential">[docs]</a><span class="k">class</span> <span class="nc">Exponential</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;Exponential solution: u&#39; = a*u + b, u(0)=A&quot;</span>

<div class="viewcode-block" id="Exponential.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">b</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">A</span><span class="o">=</span><span class="mi">1</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">a</span><span class="p">),</span> <span class="n">b</span><span class="p">,</span> <span class="n">A</span>
        <span class="k">if</span> <span class="nb">abs</span><span class="p">(</span><span class="n">a</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1E-15</span><span class="p">:</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;a=</span><span class="si">%g</span><span class="s"> is too small&#39;</span> <span class="o">%</span> <span class="n">a</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">not_suitable_solvers</span> <span class="o">=</span> <span class="p">[</span>
            <span class="s">&#39;Lsodi&#39;</span><span class="p">,</span> <span class="s">&#39;Lsodis&#39;</span><span class="p">,</span> <span class="s">&#39;Lsoibt&#39;</span><span class="p">,</span> <span class="s">&#39;MyRungeKutta&#39;</span><span class="p">,</span>
            <span class="s">&#39;MySolver&#39;</span><span class="p">,</span> <span class="s">&#39;Lsodes&#39;</span><span class="p">,</span> <span class="s">&#39;SymPy_odefun&#39;</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="Exponential.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="o">*</span><span class="n">u</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span>
</div>
<div class="viewcode-block" id="Exponential.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span>
</div>
<div class="viewcode-block" id="Exponential.f_with_args"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.f_with_args">[docs]</a>    <span class="k">def</span> <span class="nf">f_with_args</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">):</span>
        <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span> <span class="o">+</span> <span class="n">b</span>
</div>
<div class="viewcode-block" id="Exponential.f_with_kwargs"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.f_with_kwargs">[docs]</a>    <span class="k">def</span> <span class="nf">f_with_kwargs</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">b</span><span class="o">=</span><span class="mi">0</span><span class="p">):</span>
        <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span> <span class="o">+</span> <span class="n">b</span>
</div>
<div class="viewcode-block" id="Exponential.jac_with_args"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.jac_with_args">[docs]</a>    <span class="k">def</span> <span class="nf">jac_with_args</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">):</span>
        <span class="k">return</span> <span class="n">a</span>
</div>
<div class="viewcode-block" id="Exponential.jac_with_kwargs"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.jac_with_kwargs">[docs]</a>    <span class="k">def</span> <span class="nf">jac_with_kwargs</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">b</span><span class="o">=</span><span class="mi">0</span><span class="p">):</span>
        <span class="k">return</span> <span class="n">a</span>
</div>
<div class="viewcode-block" id="Exponential.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">A</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">b</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">A</span>
        <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">A</span> <span class="o">+</span> <span class="n">b</span><span class="o">/</span><span class="n">a</span><span class="p">)</span> <span class="o">-</span>  <span class="n">b</span><span class="o">/</span><span class="n">a</span>
</div>
<div class="viewcode-block" id="Exponential.default_parameters"><a class="viewcode-back" href="../../problems.html#odespy.problems.Exponential.default_parameters">[docs]</a>    <span class="k">def</span> <span class="nf">default_parameters</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="n">T</span> <span class="o">=</span> <span class="mi">5</span><span class="o">/</span><span class="nb">abs</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">)</span>
        <span class="n">d</span> <span class="o">=</span> <span class="p">{</span><span class="s">&#39;time_points&#39;</span><span class="p">:</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="mi">31</span><span class="p">),</span>
             <span class="s">&#39;atol&#39;</span><span class="p">:</span> <span class="mf">1E-4</span><span class="p">,</span> <span class="s">&#39;rtol&#39;</span><span class="p">:</span> <span class="mf">1E-3</span><span class="p">}</span>
        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">d</span><span class="p">[</span><span class="s">&#39;terminate&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">step_no</span><span class="p">:</span> <span class="n">u</span><span class="p">[</span><span class="n">step_no</span><span class="p">]</span> <span class="o">&lt;</span> <span class="bp">self</span><span class="o">.</span><span class="n">A</span><span class="o">/</span><span class="mf">20.</span>
        <span class="k">return</span> <span class="n">d</span>
</div></div>
<div class="viewcode-block" id="Logistic"><a class="viewcode-back" href="../../problems.html#odespy.problems.Logistic">[docs]</a><span class="k">class</span> <span class="nc">Logistic</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;Logistic equation&quot;</span>

<div class="viewcode-block" id="Logistic.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.Logistic.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">A</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">=</span> <span class="n">a</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">R</span> <span class="o">=</span> <span class="n">R</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="n">A</span>
</div>
<div class="viewcode-block" id="Logistic.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.Logistic.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="Logistic.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.Logistic.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">a</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u</span><span class="o">/</span><span class="n">R</span><span class="p">)</span> <span class="o">+</span> <span class="n">a</span><span class="o">*</span><span class="n">u</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="mf">1.</span><span class="o">/</span><span class="n">R</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="Logistic.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.Logistic.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">U0</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">R</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span>  <span class="c"># short form</span>
        <span class="k">return</span> <span class="n">R</span><span class="o">*</span><span class="n">U0</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">R</span> <span class="o">+</span> <span class="n">U0</span><span class="o">*</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="n">a</span><span class="o">*</span><span class="n">t</span><span class="p">)</span> <span class="o">-</span> <span class="mi">1</span><span class="p">))</span>
</div></div>
<div class="viewcode-block" id="Gaussian1"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian1">[docs]</a><span class="k">class</span> <span class="nc">Gaussian1</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;1 + Gaussian function as solution&quot;</span>

<div class="viewcode-block" id="Gaussian1.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian1.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="mf">2.0</span><span class="p">,</span> <span class="n">s</span><span class="o">=</span><span class="mf">1.0</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">s</span> <span class="o">=</span> <span class="n">c</span><span class="p">,</span> <span class="nb">float</span><span class="p">(</span><span class="n">s</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">u_exact</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="Gaussian1.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian1.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="o">-</span><span class="p">((</span><span class="n">t</span><span class="o">-</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">)</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">s</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">u</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="Gaussian1.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian1.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="o">-</span><span class="p">(</span><span class="n">t</span><span class="o">-</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">)</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">s</span><span class="o">**</span><span class="mi">2</span>
</div>
<div class="viewcode-block" id="Gaussian1.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian1.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="mi">1</span> <span class="o">+</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="mf">0.5</span><span class="o">*</span><span class="p">((</span><span class="n">t</span><span class="o">-</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">)</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">s</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
</div></div>
<div class="viewcode-block" id="Gaussian0"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian0">[docs]</a><span class="k">class</span> <span class="nc">Gaussian0</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;Gaussian function as solution&quot;</span>

<div class="viewcode-block" id="Gaussian0.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian0.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="mf">2.0</span><span class="p">,</span> <span class="n">s</span><span class="o">=</span><span class="mf">1.0</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">s</span> <span class="o">=</span> <span class="n">c</span><span class="p">,</span> <span class="nb">float</span><span class="p">(</span><span class="n">s</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">u_exact</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="Gaussian0.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian0.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="o">-</span><span class="p">((</span><span class="n">t</span><span class="o">-</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">)</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">s</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">u</span>
</div>
<div class="viewcode-block" id="Gaussian0.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian0.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="o">-</span><span class="p">(</span><span class="n">t</span><span class="o">-</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">)</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">s</span><span class="o">**</span><span class="mi">2</span>
</div>
<div class="viewcode-block" id="Gaussian0.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.Gaussian0.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="mf">0.5</span><span class="o">*</span><span class="p">((</span><span class="n">t</span><span class="o">-</span><span class="bp">self</span><span class="o">.</span><span class="n">c</span><span class="p">)</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">s</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>

</div></div>
<div class="viewcode-block" id="default_oscillator"><a class="viewcode-back" href="../../problems.html#odespy.problems.default_oscillator">[docs]</a><span class="k">def</span> <span class="nf">default_oscillator</span><span class="p">(</span><span class="n">P</span><span class="p">,</span> <span class="n">resolution_per_period</span><span class="o">=</span><span class="mi">20</span><span class="p">):</span>
    <span class="n">n</span> <span class="o">=</span> <span class="mf">4.5</span>
    <span class="n">r</span> <span class="o">=</span> <span class="n">resolution_per_period</span> <span class="c"># short form</span>
    <span class="k">print</span> <span class="s">&#39;default&#39;</span><span class="p">,</span> <span class="n">P</span><span class="p">,</span> <span class="n">n</span><span class="o">*</span><span class="n">P</span>
    <span class="n">tp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="o">*</span><span class="n">P</span><span class="p">,</span> <span class="n">r</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
    <span class="c"># atol=rtol since u approx 1, atol level set at the</span>
    <span class="c"># error RK4 produces with 20 steps per period, times 0.05</span>
    <span class="n">error_level</span> <span class="o">=</span> <span class="mf">3E-3</span>
    <span class="n">atol</span> <span class="o">=</span> <span class="n">rtol</span> <span class="o">=</span> <span class="mf">0.05</span><span class="o">*</span><span class="n">error_level</span>
    <span class="c"># E = C*dt**q</span>
    <span class="k">return</span> <span class="nb">dict</span><span class="p">(</span><span class="n">time_points</span><span class="o">=</span><span class="n">tp</span><span class="p">,</span> <span class="n">atol</span><span class="o">=</span><span class="n">atol</span><span class="p">,</span> <span class="n">rtol</span><span class="o">=</span><span class="n">rtol</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="LinearOscillator"><a class="viewcode-back" href="../../problems.html#odespy.problems.LinearOscillator">[docs]</a><span class="k">class</span> <span class="nc">LinearOscillator</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;Linear oscillator: m*u&#39;&#39; + k*u = F(t)&quot;</span>

<div class="viewcode-block" id="LinearOscillator.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.LinearOscillator.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">m</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">k</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">x0</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">v0</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">F</span><span class="o">=</span><span class="bp">None</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="p">[</span><span class="n">x0</span><span class="p">,</span> <span class="n">v0</span><span class="p">]</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">m</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">k</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">m</span><span class="p">),</span> <span class="n">k</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">F</span> <span class="o">=</span> <span class="n">F</span>
</div>
<div class="viewcode-block" id="LinearOscillator.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.LinearOscillator.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">x</span><span class="p">,</span> <span class="n">v</span> <span class="o">=</span> <span class="n">u</span>
        <span class="n">k</span><span class="p">,</span> <span class="n">m</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">k</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">m</span>
        <span class="n">F</span> <span class="o">=</span> <span class="mi">0</span> <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">F</span> <span class="ow">is</span> <span class="bp">None</span> <span class="k">else</span> <span class="bp">self</span><span class="o">.</span><span class="n">F</span><span class="p">(</span><span class="n">t</span><span class="p">)</span>
        <span class="k">return</span> <span class="p">[</span><span class="n">v</span><span class="p">,</span> <span class="o">-</span><span class="n">k</span><span class="o">/</span><span class="n">m</span><span class="o">*</span><span class="n">x</span> <span class="o">+</span> <span class="n">F</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="LinearOscillator.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.LinearOscillator.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">x</span><span class="p">,</span> <span class="n">v</span> <span class="o">=</span> <span class="n">u</span>
        <span class="n">k</span><span class="p">,</span> <span class="n">m</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">k</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">m</span>
        <span class="k">return</span> <span class="p">[[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="o">-</span><span class="n">k</span><span class="o">/</span><span class="n">m</span><span class="p">,</span> <span class="mi">0</span><span class="p">]]</span>
</div>
<div class="viewcode-block" id="LinearOscillator.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.LinearOscillator.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">F</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
            <span class="n">k</span><span class="p">,</span> <span class="n">m</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">k</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">m</span>
            <span class="n">w</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">k</span><span class="o">/</span><span class="n">m</span><span class="p">)</span>
            <span class="n">x</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">w</span><span class="o">*</span><span class="n">t</span><span class="p">)</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">w</span><span class="o">*</span><span class="n">t</span><span class="p">)</span>
            <span class="n">v</span> <span class="o">=</span> <span class="o">-</span><span class="n">w</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">U0</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">w</span><span class="o">*</span><span class="n">t</span><span class="p">)</span> <span class="o">+</span> <span class="n">w</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">U0</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">w</span><span class="o">*</span><span class="n">t</span><span class="p">)</span>
            <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">x</span><span class="p">,</span> <span class="n">v</span><span class="p">])</span><span class="o">.</span><span class="n">transpose</span><span class="p">()</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="k">return</span> <span class="bp">None</span>
</div>
<div class="viewcode-block" id="LinearOscillator.default_parameters"><a class="viewcode-back" href="../../problems.html#odespy.problems.LinearOscillator.default_parameters">[docs]</a>    <span class="k">def</span> <span class="nf">default_parameters</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">resolution_per_period</span><span class="o">=</span><span class="mi">20</span><span class="p">):</span>
        <span class="n">w</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">k</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">m</span>
        <span class="n">P</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">w</span><span class="p">)</span>  <span class="c"># period</span>
        <span class="k">return</span> <span class="n">default_oscillator</span><span class="p">(</span><span class="n">P</span><span class="p">,</span> <span class="n">resolution_per_period</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="LinearOscillator.kinetic_energy"><a class="viewcode-back" href="../../problems.html#odespy.problems.LinearOscillator.kinetic_energy">[docs]</a>    <span class="k">def</span> <span class="nf">kinetic_energy</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">):</span>
        <span class="n">v</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span>
        <span class="k">return</span> <span class="mf">0.5</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">m</span><span class="o">*</span><span class="n">v</span><span class="o">**</span><span class="mi">2</span>
</div>
<div class="viewcode-block" id="LinearOscillator.potential_energy"><a class="viewcode-back" href="../../problems.html#odespy.problems.LinearOscillator.potential_energy">[docs]</a>    <span class="k">def</span> <span class="nf">potential_energy</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">):</span>
        <span class="n">x</span> <span class="o">=</span> <span class="n">u</span><span class="p">[:,</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">return</span> <span class="mf">0.5</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">k</span><span class="o">*</span><span class="n">x</span><span class="o">**</span><span class="mi">2</span>

<span class="c"># offer analytical solution with damping, also with sin/cos excitation</span>
</div></div>
<div class="viewcode-block" id="ComplexOscillator"><a class="viewcode-back" href="../../problems.html#odespy.problems.ComplexOscillator">[docs]</a><span class="k">class</span> <span class="nc">ComplexOscillator</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Harmonic oscillator (u&#39;&#39; + w*u = 0) expressed as a complex</span>
<span class="sd">    ODE: u&#39; = i*w*u.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;Complex oscillator: u&#39; = i*w*u&quot;</span>

<div class="viewcode-block" id="ComplexOscillator.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.ComplexOscillator.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">w</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">U0</span><span class="o">=</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">]):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">w</span> <span class="o">=</span> <span class="n">w</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">not_suitable_solvers</span> <span class="o">=</span> <span class="n">list_not_suitable_complex_solvers</span><span class="p">()</span>
</div>
<div class="viewcode-block" id="ComplexOscillator.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.ComplexOscillator.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="mi">1</span><span class="n">j</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">w</span><span class="o">*</span><span class="n">u</span>
</div>
<div class="viewcode-block" id="ComplexOscillator.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.ComplexOscillator.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="mi">1</span><span class="n">j</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">w</span>
</div>
<div class="viewcode-block" id="ComplexOscillator.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.ComplexOscillator.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="mi">1</span><span class="n">j</span><span class="o">*</span><span class="bp">self</span><span class="o">.</span><span class="n">w</span><span class="o">*</span><span class="n">t</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="ComplexOscillator.default_parameters"><a class="viewcode-back" href="../../problems.html#odespy.problems.ComplexOscillator.default_parameters">[docs]</a>    <span class="k">def</span> <span class="nf">default_parameters</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">resolution_per_period</span><span class="o">=</span><span class="mi">20</span><span class="p">):</span>
        <span class="n">P</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">w</span><span class="p">)</span>  <span class="c"># period</span>
        <span class="k">return</span> <span class="n">default_oscillator</span><span class="p">(</span><span class="n">P</span><span class="p">,</span> <span class="n">resolution_per_period</span><span class="p">)</span>
</div></div>
<div class="viewcode-block" id="StiffSystem2x2"><a class="viewcode-back" href="../../problems.html#odespy.problems.StiffSystem2x2">[docs]</a><span class="k">class</span> <span class="nc">StiffSystem2x2</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;Potentially stiff 2x2 system u&#39; = 1/eps*u&quot;</span>
    <span class="n">stiff</span> <span class="o">=</span> <span class="bp">True</span>

<div class="viewcode-block" id="StiffSystem2x2.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.StiffSystem2x2.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">eps</span><span class="o">=</span><span class="mf">1E-2</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">eps</span> <span class="o">=</span> <span class="n">eps</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>
        <span class="k">if</span> <span class="mf">0.2</span> <span class="o">&lt;=</span> <span class="n">eps</span> <span class="o">&lt;=</span> <span class="mi">5</span><span class="p">:</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">stiff</span> <span class="o">=</span> <span class="bp">False</span>
</div>
<div class="viewcode-block" id="StiffSystem2x2.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.StiffSystem2x2.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="p">[</span><span class="o">-</span><span class="n">u</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">eps</span><span class="o">*</span><span class="n">u</span><span class="p">[</span><span class="mi">1</span><span class="p">]]</span>
</div>
<div class="viewcode-block" id="StiffSystem2x2.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.StiffSystem2x2.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="p">[[</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">eps</span><span class="p">]]</span>
</div>
<div class="viewcode-block" id="StiffSystem2x2.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.StiffSystem2x2.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">t</span><span class="p">),</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="mf">1.</span><span class="o">/</span><span class="bp">self</span><span class="o">.</span><span class="n">eps</span><span class="o">*</span><span class="n">t</span><span class="p">)])</span><span class="o">.</span><span class="n">transpose</span><span class="p">()</span>

</div></div>
<div class="viewcode-block" id="VanDerPolOscillator"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator">[docs]</a><span class="k">class</span> <span class="nc">VanDerPolOscillator</span><span class="p">(</span><span class="n">Problem</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Classical van der Pool oscillator:</span>

<span class="sd">    .. math:: y&#39;&#39; = \mu (1 - y^2) y&#39; - y</span>

<span class="sd">    with initial conditions :math:`y(0)=2, y&#39;(0)=1`.</span>
<span class="sd">    The equation is rewritten as a system</span>

<span class="sd">    .. math::</span>
<span class="sd">             u_0&#39; &amp;= u_1, \\</span>
<span class="sd">             u_1&#39; &amp;= \mu (1-u_0^2)u_1 - u_0</span>

<span class="sd">    with a Jacobian</span>

<span class="sd">    .. math::</span>
<span class="sd">             \left(\begin{array}{cc}</span>
<span class="sd">             0 &amp; 1\\</span>
<span class="sd">             -2\mu u_0 - 1 &amp; \mu (1-u_0^2)</span>
<span class="sd">             \end{array}\right)</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">short_description</span> <span class="o">=</span> <span class="s">&quot;Van der Pol oscillator&quot;</span>

<div class="viewcode-block" id="VanDerPolOscillator.__init__"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.__init__">[docs]</a>    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">U0</span><span class="o">=</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="n">mu</span><span class="o">=</span><span class="mf">3.</span><span class="p">,</span> <span class="n">f77</span><span class="o">=</span><span class="bp">False</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">U0</span> <span class="o">=</span> <span class="n">U0</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">mu</span> <span class="o">=</span> <span class="n">mu</span>

        <span class="c"># Compile F77 functions</span>
        <span class="k">if</span> <span class="n">f77</span><span class="p">:</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">f_f77</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">jac_f77_radau5</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">jac_f77_lsode</span> <span class="o">=</span> \
                        <span class="n">compile_f77</span><span class="p">([</span><span class="bp">self</span><span class="o">.</span><span class="n">str_f_f77</span><span class="p">(),</span>
                                     <span class="bp">self</span><span class="o">.</span><span class="n">str_jac_f77_radau5</span><span class="p">(),</span>
                                     <span class="bp">self</span><span class="o">.</span><span class="n">str_jac_f77_lsode</span><span class="p">()])</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.f"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.f">[docs]</a>    <span class="k">def</span> <span class="nf">f</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">u_0</span><span class="p">,</span> <span class="n">u_1</span> <span class="o">=</span> <span class="n">u</span>
        <span class="n">mu</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">mu</span>
        <span class="k">return</span> <span class="p">[</span><span class="n">u_1</span><span class="p">,</span> <span class="n">mu</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u_0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">u_1</span> <span class="o">-</span> <span class="n">u_0</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.jac"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.jac">[docs]</a>    <span class="k">def</span> <span class="nf">jac</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="n">u_0</span><span class="p">,</span> <span class="n">u_1</span> <span class="o">=</span> <span class="n">u</span>
        <span class="n">mu</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">mu</span>
        <span class="k">return</span> <span class="p">[[</span><span class="mf">0.</span><span class="p">,</span> <span class="mf">1.</span><span class="p">],</span>
                <span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">mu</span><span class="o">*</span><span class="n">u_0</span><span class="o">*</span><span class="n">u_1</span> <span class="o">-</span> <span class="mi">1</span><span class="p">,</span> <span class="n">mu</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u_0</span><span class="o">**</span><span class="mi">2</span><span class="p">)]]</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.f_with_args"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.f_with_args">[docs]</a>    <span class="k">def</span> <span class="nf">f_with_args</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">mu</span><span class="p">):</span>
        <span class="n">u_0</span><span class="p">,</span> <span class="n">u_1</span> <span class="o">=</span> <span class="n">u</span>
        <span class="k">return</span> <span class="p">[</span><span class="n">u_1</span><span class="p">,</span> <span class="n">mu</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u_0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">u_1</span> <span class="o">-</span> <span class="n">u_0</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.jac_with_args"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.jac_with_args">[docs]</a>    <span class="k">def</span> <span class="nf">jac_with_args</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">mu</span><span class="p">):</span>
        <span class="n">u_0</span><span class="p">,</span> <span class="n">u_1</span> <span class="o">=</span> <span class="n">u</span>
        <span class="k">return</span> <span class="p">[[</span><span class="mf">0.</span><span class="p">,</span> <span class="mf">1.</span><span class="p">],</span>
                <span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">mu</span><span class="o">*</span><span class="n">u_0</span><span class="o">*</span><span class="n">u_1</span> <span class="o">-</span> <span class="mi">1</span><span class="p">,</span> <span class="n">mu</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u_0</span><span class="o">**</span><span class="mi">2</span><span class="p">)]]</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.f_with_kwargs"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.f_with_kwargs">[docs]</a>    <span class="k">def</span> <span class="nf">f_with_kwargs</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">mu</span><span class="o">=</span><span class="mi">3</span><span class="p">):</span>
        <span class="n">u_0</span><span class="p">,</span> <span class="n">u_1</span> <span class="o">=</span> <span class="n">u</span>
        <span class="k">return</span> <span class="p">[</span><span class="n">u_1</span><span class="p">,</span> <span class="n">mu</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u_0</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">u_1</span> <span class="o">-</span> <span class="n">u_0</span><span class="p">]</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.jac_with_kwargs"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.jac_with_kwargs">[docs]</a>    <span class="k">def</span> <span class="nf">jac_with_kwargs</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">mu</span><span class="o">=</span><span class="mi">3</span><span class="p">):</span>
        <span class="n">u_0</span><span class="p">,</span> <span class="n">u_1</span> <span class="o">=</span> <span class="n">u</span>
        <span class="k">return</span> <span class="p">[[</span><span class="mf">0.</span><span class="p">,</span> <span class="mf">1.</span><span class="p">],</span>
                <span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">mu</span><span class="o">*</span><span class="n">u_0</span><span class="o">*</span><span class="n">u_1</span> <span class="o">-</span> <span class="mi">1</span><span class="p">,</span> <span class="n">mu</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">u_0</span><span class="o">**</span><span class="mi">2</span><span class="p">)]]</span>

</div>
<div class="viewcode-block" id="VanDerPolOscillator.str_f_f77"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.str_f_f77">[docs]</a>    <span class="k">def</span> <span class="nf">str_f_f77</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Return f(u,t) as Fortran source code string.&quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="s">&quot;&quot;&quot;</span>
<span class="s">      subroutine f_f77(neq, t, u, udot)</span>
<span class="s">Cf2py intent(hide) neq</span>
<span class="s">Cf2py intent(out) udot</span>
<span class="s">      integer neq</span>
<span class="s">      double precision t, u, udot</span>
<span class="s">      dimension u(neq), udot(neq)</span>
<span class="s">      udot(1) = u(2)</span>
<span class="s">      udot(2) = </span><span class="si">%g</span><span class="s">*(1 - u(1)**2)*u(2) - u(1)</span>
<span class="s">      return</span>
<span class="s">      end</span>
<span class="s">&quot;&quot;&quot;</span> <span class="o">%</span> <span class="bp">self</span><span class="o">.</span><span class="n">mu</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.str_jac_f77_fadau5"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.str_jac_f77_fadau5">[docs]</a>    <span class="k">def</span> <span class="nf">str_jac_f77_fadau5</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Return Jacobian Fortran source code string for Radau5.&quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="s">&quot;&quot;&quot;</span>
<span class="s">      subroutine jac_f77_radau5(neq,t,u,dfu,ldfu,rpar,ipar)</span>
<span class="s">Cf2py intent(hide) neq,rpar,ipar</span>
<span class="s">Cf2py intent(in)   t,u,ldfu</span>
<span class="s">Cf2py intent(out) dfu</span>
<span class="s">      integer neq,ipar,ldfu</span>
<span class="s">      double precision t,u,dfu,rpar</span>
<span class="s">      dimension u(neq),dfu(ldfu,neq),rpar(*),ipar(*)</span>
<span class="s">      dfu(1,1) = 0</span>
<span class="s">      dfu(1,2) = 1</span>
<span class="s">      dfu(2,1) = -2*</span><span class="si">%g</span><span class="s">*u(1)*u(2) - 1</span>
<span class="s">      dfu(2,2) = </span><span class="si">%g</span><span class="s">*(1-u(1)**2)</span>
<span class="s">      return</span>
<span class="s">      end</span>
<span class="s">&quot;&quot;&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">mu</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">mu</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.str_jac_f77_lsode_dense"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.str_jac_f77_lsode_dense">[docs]</a>    <span class="k">def</span> <span class="nf">str_jac_f77_lsode_dense</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Return Fortran source for dense Jacobian matrix in LSODE format.&quot;&quot;&quot;</span>
        <span class="k">return</span> <span class="s">&quot;&quot;&quot;</span>
<span class="s">      subroutine jac_f77(neq, t, u, ml, mu, pd, nrowpd)</span>
<span class="s">Cf2py intent(hide) neq, ml, mu, nrowpd</span>
<span class="s">Cf2py intent(out) pd</span>
<span class="s">      integer neq, ml, mu, nrowpd</span>
<span class="s">      double precision t, u, pd</span>
<span class="s">      dimension u(neq), pd(nrowpd,neq)</span>
<span class="s">      pd(1,1) = 0</span>
<span class="s">      pd(1,2) = 1</span>
<span class="s">      pd(2,1) = -2*</span><span class="si">%g</span><span class="s">*u(1)*u(2) - 1</span>
<span class="s">      pd(2,2) = </span><span class="si">%g</span><span class="s">*(1 - u(1)**2)</span>
<span class="s">      return</span>
<span class="s">      end</span>
<span class="s">&quot;&quot;&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">mu</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">mu</span><span class="p">)</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.u_exact"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.u_exact">[docs]</a>    <span class="k">def</span> <span class="nf">u_exact</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">t</span><span class="p">):</span>
        <span class="k">if</span> <span class="nb">abs</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">mu</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1E-14</span><span class="p">:</span>
            <span class="n">x</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">t</span><span class="p">)</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">t</span><span class="p">)</span>
            <span class="n">v</span> <span class="o">=</span> <span class="o">-</span><span class="bp">self</span><span class="o">.</span><span class="n">U0</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">t</span><span class="p">)</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">U0</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">t</span><span class="p">)</span>
            <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">x</span><span class="p">,</span> <span class="n">v</span><span class="p">])</span><span class="o">.</span><span class="n">transpose</span><span class="p">()</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="k">return</span> <span class="bp">None</span>
</div>
<div class="viewcode-block" id="VanDerPolOscillator.default_parameters"><a class="viewcode-back" href="../../problems.html#odespy.problems.VanDerPolOscillator.default_parameters">[docs]</a>    <span class="k">def</span> <span class="nf">default_parameters</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">resolution_per_period</span><span class="o">=</span><span class="mi">20</span><span class="p">):</span>
        <span class="c"># Period = 2*pi for mu=0</span>
        <span class="n">P</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">pi</span>
        <span class="k">return</span> <span class="n">default_oscillator</span><span class="p">(</span><span class="n">P</span><span class="p">,</span> <span class="n">resolution_per_period</span><span class="p">)</span>


</div></div>
<div class="viewcode-block" id="tester"><a class="viewcode-back" href="../../problems.html#odespy.problems.tester">[docs]</a><span class="k">def</span> <span class="nf">tester</span><span class="p">(</span><span class="n">problems</span><span class="p">,</span> <span class="n">methods</span><span class="p">,</span> <span class="n">time_points</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span> <span class="n">compare_tol</span><span class="o">=</span><span class="mf">1E-4</span><span class="p">,</span>
           <span class="n">solver_prm</span><span class="o">=</span><span class="p">{}):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    `problems` is a list of Problem subclass instances made ready.</span>
<span class="sd">    `methods` is a list of strings holding the method names.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="kn">import</span> <span class="nn">odespy</span>
    <span class="n">results</span> <span class="o">=</span> <span class="p">{}</span>
    <span class="n">error_msg</span> <span class="o">=</span> <span class="p">{}</span>
    <span class="k">for</span> <span class="n">problem</span> <span class="ow">in</span> <span class="n">problems</span><span class="p">:</span>
        <span class="n">pname</span> <span class="o">=</span> <span class="n">problem</span><span class="o">.</span><span class="n">__class__</span><span class="o">.</span><span class="n">__name__</span>
        <span class="k">print</span> <span class="s">&#39;problem &#39;</span><span class="p">,</span> <span class="n">pname</span>
        <span class="n">methods4problem</span> <span class="o">=</span> <span class="p">[</span><span class="n">method</span> <span class="k">for</span> <span class="n">method</span> <span class="ow">in</span> <span class="n">methods</span>
                           <span class="k">if</span> <span class="n">method</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">problem</span><span class="o">.</span><span class="n">not_suitable_solvers</span><span class="p">]</span>
        <span class="n">defaults</span> <span class="o">=</span> <span class="n">problem</span><span class="o">.</span><span class="n">default_parameters</span><span class="p">()</span>
        <span class="k">if</span> <span class="n">time_points</span> <span class="ow">is</span> <span class="bp">None</span><span class="p">:</span>
            <span class="n">time_points</span> <span class="o">=</span> <span class="n">defaults</span><span class="p">[</span><span class="s">&#39;time_points&#39;</span><span class="p">]</span>
        <span class="n">results</span><span class="p">[</span><span class="n">pname</span><span class="p">]</span> <span class="o">=</span> <span class="p">{</span><span class="s">&#39;t&#39;</span><span class="p">:</span> <span class="n">time_points</span><span class="p">}</span>
        <span class="n">error_msg</span><span class="p">[</span><span class="n">pname</span><span class="p">]</span> <span class="o">=</span> <span class="p">{}</span>

        <span class="n">this_solver_prm</span> <span class="o">=</span> <span class="n">solver_prm</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
        <span class="n">names_in_problem</span> <span class="o">=</span> <span class="s">&#39;jac&#39;</span><span class="p">,</span> <span class="s">&#39;spcrad&#39;</span><span class="p">,</span> <span class="s">&#39;u_exact&#39;</span><span class="p">,</span>
        <span class="k">for</span> <span class="n">name</span> <span class="ow">in</span> <span class="n">names_in_problem</span><span class="p">:</span>
            <span class="c"># Set parameter if problem has it and user has</span>
            <span class="c"># not specified it in solver_prm</span>
            <span class="k">if</span> <span class="n">name</span> <span class="ow">in</span> <span class="n">problem</span> <span class="ow">and</span> <span class="n">name</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">solver_prm</span><span class="p">:</span>
                <span class="n">this_solver_prm</span><span class="p">[</span><span class="n">name</span><span class="p">]</span> <span class="o">=</span> <span class="nb">getattr</span><span class="p">(</span><span class="n">problem</span><span class="p">,</span> <span class="n">name</span><span class="p">)</span>
        <span class="n">names_in_defaults</span> <span class="o">=</span> <span class="s">&#39;atol&#39;</span><span class="p">,</span> <span class="s">&#39;rtol&#39;</span><span class="p">,</span> <span class="s">&#39;max_iter&#39;</span><span class="p">,</span> <span class="s">&#39;eps_iter&#39;</span>
        <span class="k">for</span> <span class="n">name</span> <span class="ow">in</span> <span class="n">names_in_defaults</span><span class="p">:</span>
            <span class="k">if</span> <span class="n">name</span> <span class="ow">in</span> <span class="n">defaults</span> <span class="ow">and</span> <span class="n">name</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">solver_prm</span><span class="p">:</span>
                <span class="n">this_solver_prm</span><span class="p">[</span><span class="n">name</span><span class="p">]</span> <span class="o">=</span> <span class="n">defaults</span><span class="p">[</span><span class="n">name</span><span class="p">]</span>

        <span class="k">for</span> <span class="n">method</span> <span class="ow">in</span> <span class="n">methods4problem</span><span class="p">:</span>
            <span class="k">print</span> <span class="s">&#39;  testing&#39;</span><span class="p">,</span> <span class="n">method</span><span class="p">,</span>
            <span class="n">solver</span> <span class="o">=</span> <span class="nb">eval</span><span class="p">(</span><span class="s">&#39;odespy.&#39;</span><span class="o">+</span><span class="n">method</span><span class="p">)(</span><span class="n">problem</span><span class="o">.</span><span class="n">f</span><span class="p">)</span>

            <span class="c"># Important to set parameters before setting initial cond.</span>
            <span class="n">solver</span><span class="o">.</span><span class="n">set</span><span class="p">(</span><span class="o">**</span><span class="n">this_solver_prm</span><span class="p">)</span>

            <span class="n">solver</span><span class="o">.</span><span class="n">set_initial_condition</span><span class="p">(</span><span class="n">problem</span><span class="o">.</span><span class="n">get_initial_condition</span><span class="p">())</span>

            <span class="n">u</span><span class="p">,</span> <span class="n">t</span> <span class="o">=</span> <span class="n">solver</span><span class="o">.</span><span class="n">solve</span><span class="p">(</span><span class="n">time_points</span><span class="p">,</span> <span class="n">problem</span><span class="o">.</span><span class="n">terminate</span><span class="p">)</span>

            <span class="n">error</span> <span class="o">=</span> <span class="n">problem</span><span class="o">.</span><span class="n">verify</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">t</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">error</span> <span class="ow">is</span> <span class="ow">not</span> <span class="bp">None</span><span class="p">:</span>
                <span class="n">results</span><span class="p">[</span><span class="n">pname</span><span class="p">][</span><span class="n">method</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">error</span><span class="p">)</span>
                <span class="k">print</span> <span class="n">error</span><span class="p">,</span>
                <span class="k">if</span> <span class="n">error</span> <span class="o">&gt;</span> <span class="n">compare_tol</span><span class="p">:</span>
                    <span class="k">print</span> <span class="s">&#39;WARNING: tolerance </span><span class="si">%.0E</span><span class="s"> exceeded&#39;</span> <span class="o">%</span> <span class="n">compare_tol</span><span class="p">,</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">results</span><span class="p">[</span><span class="n">pname</span><span class="p">][</span><span class="n">method</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">u</span><span class="p">,)</span>
            <span class="k">print</span>
    <span class="k">return</span> <span class="n">results</span>
</div>
<span class="sd">&quot;&quot;&quot;</span>

<span class="sd">Radau5 {&#39;rtol&#39;: 0.01, &#39;atol&#39;: 0.01, &#39;max_iter&#39;: 10, &#39;jac&#39;: &lt;bound method Linear2.jac of &lt;problems.Linear2 instance at 0x1fdd8c0&gt;&gt;, &#39;eps_iter&#39;: 1e-06, &#39;spcrad&#39;: &lt;bound method Linear2.spcrad of &lt;problems.Linear2 instance at 0x1fdd8c0&gt;&gt;}</span>
<span class="sd">1.7763568394e-15</span>
<span class="sd">Radau5Explicit {&#39;rtol&#39;: 0.01, &#39;atol&#39;: 0.01, &#39;max_iter&#39;: 10, &#39;jac&#39;: &lt;bound method Linear2.jac of &lt;problems.Linear2 instance at 0x1fdd8c0&gt;&gt;, &#39;eps_iter&#39;: 1e-06, &#39;spcrad&#39;: &lt;bound method Linear2.spcrad of &lt;problems.Linear2 instance at 0x1fdd8c0&gt;&gt;}</span>
<span class="sd">1.7763568394e-15</span>
<span class="sd">Radau5Implicit {&#39;rtol&#39;: 0.01, &#39;atol&#39;: 0.01, &#39;max_iter&#39;: 10, &#39;jac&#39;: &lt;bound method Linear2.jac of &lt;problems.Linear2 instance at 0x1fdd8c0&gt;&gt;, &#39;eps_iter&#39;: 1e-06, &#39;spcrad&#39;: &lt;bound method Linear2.spcrad of &lt;problems.Linear2 instance at 0x1fdd8c0&gt;&gt;}</span>

<span class="sd">    # Compare with exact solution</span>
<span class="sd">    # Compare solutions within tolerance</span>
<span class="sd">    failure = {}</span>
<span class="sd">    for pname in problems:</span>
<span class="sd">        failure[pname] = {}</span>
<span class="sd">        for i, mname in enumeate(methods):</span>
<span class="sd">            if i == 0:</span>
<span class="sd">                reference_solution, t = results[(mname, pname)]</span>
<span class="sd">                reference_name = mname</span>
<span class="sd">            else:</span>
<span class="sd">                u = results[(mname, pname)][0]</span>
<span class="sd">                r = np.allclose(reference_solution, u,</span>
<span class="sd">                                compare_tol, 0.1*compare_tol)</span>
<span class="sd">                if r:</span>
<span class="sd">                    failure[pname][mname] = False</span>
<span class="sd">                else:</span>
<span class="sd">                    failure[pname][mname] = np.abs(reference_solution - u).max()</span>
<span class="sd">    return failure</span>
<span class="sd">&quot;&quot;&quot;</span>
</pre></div>

          </div>
        </div>
      </div>
      <div class="sphinxsidebar">
        <div class="sphinxsidebarwrapper">
<div id="searchbox" style="display: none">
  <h3>Quick search</h3>
    <form class="search" action="../../search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="../../genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="../../py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="../../np-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li><a href="../../index.html">Odespy API 0.2 documentation</a> &raquo;</li>
          <li><a href="../index.html" >Module code</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer">
        &copy; Copyright 2012, Liwei Wang and Hans Petter Langtangen.
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.1.2.
    </div>
  </body>
</html>